clear all
use "datasets/MRA_waste_sites"

*-------------------------------------------------------------------------*
* FIGURE A3: Histogram of effect size.
*-------------------------------------------------------------------------*

histogram elas if elas>=-1 & elas<=1, xtitle(Effect size) scheme(lean1) fraction width (0.1)
graph export "figures/Figure A3. Histogram of effect size.tif", replace as(tif) ///4 observations are outside the frame of this plot

*-------------------------------------------------------------------------*
* TABLE 1: Summary statistics of moderators.
*-------------------------------------------------------------------------*

*Continuous and dummy variables
estpost tabstat elas_SE publish year_publish site_m HDI GDP_pc_national_2010_USD NPL_n NPL_y not_USA NPL_u active_n active_y active_u job_n job_y job_u haz non_haz nuclear soil air water element_u North_America Asia Europe Other_continent cleanup0 cleanup1 cleanup2 cleanup3 cleanup_u dist_mean miles_km dist_greater_mean sample sale_ind num_sig_var num_expl num_struc num_nb num_env oth_disamen oth_amen access industry demoecon time_control direction interaction log_log OLS_spatial, quietly columns(statistics) statistics (mean sd count)

esttab using "tables/summary moderators", replace cells("mean(fmt(2)) sd(fmt(2)) count(fmt(0))") noobs nonumbers rtf addnotes("Notes: The asterisk denotes moderators that have not been used in this field before.")

* Number of observations for the levels of categorical variables
tab1 NPL active job site_cat element region cleanup_stage,sort

*-------------------------------------------------------------------------*
* TABLE 2: Distribution of effect size by selected study characteristics.
*-------------------------------------------------------------------------*

*Additional information on the number of studies that exhibit a particular feature
egen tag_study = tag(ID_Study) // leads to wrong numbers of studies per category if the variable varies within the study e.g. clean-up status / site_cat --> manual counting 

tabstat elas, statistics (n mean p5 p95)
tabstat elas if tag_study, statistics (n)

bys site_cat: tabstat elas, statistics (n mean p5 p95)
bys site_cat: tabstat elas if tag_study, statistics (n)

bys element: tabstat elas , statistics (n mean p5 p95)
bys element: tabstat elas if tag_study, statistics (n)

bys region : tabstat elas, statistics (n mean p5 p95)
bys region : tabstat elas if tag_study, statistics (n)
bys Other_continent: tabstat elas, statistics (n mean p5 p95)

bys cleanup_stage: tabstat elas, statistics (n mean p5 p95)
bys cleanup_stage: tabstat elas if tag_study, statistics (n)

bys publish: tabstat elas, statistics (n mean p5 p95)
bys publish: tabstat elas if tag_study, statistics (n)

*--------------------------------------------------------------------------------*
* FIGURE A4: Funnel Plot for individual estimates and study means
*--------------------------------------------------------------------------------*
*by estimate
metafunnel elas elas_SE, by(site_cat) scheme(s2mono) graphregion(color(white))subtitle("Funnel plot: All estimates by site catagory") ylab(1.5(0.5)0, nogrid) xtitle(Effect size) ytitle(Standard error)
graph rename Graph funnel_all_estimates, replace

*Creating the study means
by ID_Study, sort: egen meanelas = mean(elas)
by ID_Study, sort: egen meanelas_SE = mean(elas_SE)
bys ID_Study: gen one_per_study= 1 if _n == 1 

*by study means
metafunnel meanelas meanelas_SE if one_per_study==1, by(site_cat) scheme(s2mono) graphregion(color(white)) xtitle(Effect size) ytitle(Standard error) ylab(0.3(0.1)0, nogrid) subtitle("Funnel plot: Study means by site catagory")  
graph rename Graph funnel_mean_estimates, replace

graph combine funnel_all_estimates funnel_mean_estimates ,graphregion(color(white))
graph export "figures/Figure A4. Funnel plots.tif", replace as(tif)

*--------------------------------------------------------------------------------*
* TABLE 3: Publication bias                                                          *
*--------------------------------------------------------------------------------*

*RANDOM EFFECTS

*No publication bias control: WLS-RE
quietly metareg elas, wsse(elas_SE) tau2test  // Q=13511 p=0.000
scalar tau2 = e(tau2)
display tau2
gen elas_var_re= elas_var + tau2
gen elas_SE_re = sqrt(elas_var_re)
gen precision_re= 1/elas_SE_re

foreach var in elas elas_SE elas_var{
gen `var'1 = `var'/elas_SE_re
}
reg elas1 precision_re, nocons cluster(ID_Study)
outreg2 using "tables/Table 3. Publication bias and genuine effect size.rtf", replace ctitle(WLS-RE) nodepvar stats(coef se) dec(3) 

* Linear publication bias control: WLS-RE FAT-PET
quietly metareg elas elas_SE, wsse(elas_SE) tau2test  // Q=11350 p=0.000
scalar tau2 = e(tau2)
display tau2
gen elas_var_re_pb= elas_var + tau2
gen elas_SE_re_pb = sqrt(elas_var_re_pb)
gen precision_re_pb= 1/(elas_SE_re_pb)

foreach var in elas elas_SE elas_var{
gen `var'2 = `var'/elas_SE_re_pb
}
reg elas2 precision_re_pb elas_SE2, nocons cluster(ID_Study)
outreg2 using "tables/Table 3. Publication bias and genuine effect size.rtf", append ctitle(WLS-RE FAT-PET) nodepvar stats(coef se) dec(3)  

*Squared publication bias control: WLS-RE PEESE
quietly metareg elas elas_var, wsse(elas_SE) tau2test  // Q=13448 p=0.000
scalar tau2 = e(tau2)
display tau2
gen elas_var_re_pb_sq= elas_var + tau2
gen elas_SE_re_pb_sq = sqrt(elas_var_re_pb_sq)
gen precision_re_pb_sq= 1/(elas_SE_re_pb_sq)

foreach var in elas elas_SE elas_var{
gen `var'3 = `var'/elas_SE_re_pb_sq
}
reg elas3 precision_re_pb_sq elas_var3, nocons cluster(ID_Study)
outreg2 using "tables/Table 3. Publication bias and genuine effect size.rtf", append ctitle(WLS-RE PEESE) nodepvar stats(coef se) dec(3) 

*FIXED EFFECTS

*generate the necessary FE weighted observations
foreach var in elas elas_SE elas_var{
gen `var'0 = `var'*log_precision
}

// no publication bias control: WLS-FE
reg elas0 log_precision, nocons cluster(ID_Study)
outreg2 using "tables/Table 3. Publication bias and genuine effect size.rtf", append ctitle(WLS-FE) nodepvar stats(coef se) dec(3) 

// Linear publication bias control: WLS-FE FAT-PET
reg elas0 log_precision elas_SE0, nocons cluster(ID_Study)
outreg2 using "tables/Table 3. Publication bias and genuine effect size.rtf", append ctitle(WLS-FE FAT-PET) nodepvar stats(coef se) dec(3) 

// Squared publication bias control: WLS-fE PEESE
reg elas0 log_precision elas_var0, nocons cluster(ID_Study)
outreg2 using "tables/Table 3. Publication bias and genuine effect size.rtf", append ctitle(WLS-FE PEESE) nodepvar stats(coef se) dec(3) 

*------------------------------------------------------------------------*
*Table B1. Results for alternative tests of publication selection bias
*------------------------------------------------------------------------*

*Weighted average of the adequately powered (WAAP)

reg elas0 log_precision, nocons // This is the WLS-FE simple meta-average used as starting point
display _se[log_precision]
gen threshold = _se[log_precision]/2.8, // this generates the threshold for adequate power (80%) with a significance level of 0.05
reg elas0 log_precision if elas_SE <= threshold, nocons // This is the WAAP 

*Stanley (2010) 'Top 10'

sum elas_SE, detail
local ten_percent=r(p10)
display `ten_percent'
reg elas if elas_SE <= `ten_percent'  

* Results for Andrews & Kasy (2019) and Furukawa (2019) were calculated using the online-tools respectively, see
* https://maxkasy.github.io/home/metastudy and https://github.com/Chishio318/stem-based_method

clear

